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Abstract 

We consider some fractional extensions of the recursive differential equation governing the 
Poisson process, i.e. 

^P*(t) = -A(p fc (t) -p fc -i(0), k > l,t > 

by introducing fractional time-derivatives of order z/, 2i/,...,nv . We show that the so-called 
"Generalized Mittag-Leffler functions" p(x) (introduced by Prabhakar [20]) arise as solu- 
tions of these equations. The corresponding processes are proved to be renewal, with density 
of the intearrival times (represented by Mittag-Leffler functions) possessing power, instead of 
exponential, decay, for t — > oo. On the other hand, near the origin the behavior of the law of 
the interarrival times drastically changes for the parameter v varying in (0, 1] . 

For integer values of v, these models can be viewed as a higher-order Poisson processes, 
connected with the standard case by simple and explict relationships. 

Key words: Fractional difference-differential equations; Generalized Mittag-Leffler func- 
tions; Fractional Poisson processes; Processes with random time; Renewal function; Cox pro- 
cess. 
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1 Introduction 

Many well-known differential equations have been extended by introducing fractional-order deriva- 
tives with respect to time (for instance, the heat, wave and telegraph equations, as well as the 
higher-order heat-type equations) or with respect to space (for instance, the equations involving the 
Riesz fractional operator). 

Fractional versions of the Poisson processes have been already presented and studied in the 
literature: in [9] the so-called fractional master equation was considered. A similar model was treated 
in [12], where the equation governing the probability distribution of the homogeneous Poisson process 
was modified, by introducing the Riemann-Liouville fractional derivative. The results are given in 
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analytical form, in terms of infinite series or successive derivatives of Mittag-Leffler functions. We 
recall the definition of the (two-parameters) Mittag-Leffler function: 



E a ^(x) = y2 , a,f3eC, Re(a),Re(/3)>0,xeR, (1.1) 

r— 

(see [20], 1.2). 

A different definition of Poisson fractional process has been proposed by Wang and Wen [30] and 
successively studied in [31]-[32]: in analogy to the well-known fractional Brownian motion, the new 
process is defined as a stochastic integral with respect to the Poisson measure. It displays properties 
similar to fractional Brownian motion, such as self-similarity (in the wide-sense) and long-range 
dependence. 

Another approach was followed by Repin and Saichev [23]: they start by generalizing, in a 
fractional sense, the distribution of the interarrival times Uj between two Poisson events. This is 
expressed, in terms of Mittag-Leffler functions, as follows, for v 6 (0,1] : 

f(t) =Pr{U j e dt}/dt=-jE vA (-n = { L\ , *>0 (1.2) 

m=l ^ ' 

and coincides with the solution to the fractional equation 

d v f(t) 



dt v 



= -f(t) + S(t), t>0 (1.3) 



where <$(•) denotes the Dirac delta function and again the fractional derivative is intended in the 
Riemann-Liouville sense. For v — 1 formula ( |1.2[ ) reduces to the well-known density appearing in 
the case of a homogeneous Poisson process, N(t),t > 0, with intensity A = 1, i.e. f(t) — e~*. 

The same approach is followed by Mainardi et al. [14]- [15]- [16], where a deep analysis of the 
related process is performed: it turns out to be a true renewal process, while it looses the Markovian 



property. Their first step is the study of the following fractional equation (instead of ( 1 . 3 1 ) 

^ = -*«>, M 

with initial condition ip(® + ) — 1 an d with fractional derivative defined in the Caputo sense. The 



solution ip(t) — E u _i(—t v ) to (1.4) represents the survival probability of the fractional Poisson 
process. As a consequence its probability distribution is expressed in terms of derivatives of Mittag- 
Leffler functions, while the density of the fc-th event waiting time is a fractional generalization of 
the Erlang distribution and coincides with the fc-fold convolution of \\.2\ . 

The analysis carried out by Beghin and Orsingher [2] starts, as in [11], from the generalization of 
the equation governing the Poisson process, where the time-derivative is substituted by the fractional 
derivative (in the Caputo sense) of order v £ (0,1]: 

^. = -\(p k -p k _ 1 ), fc>0, (1.5) 

with initial conditions 

»(°)= k>l 
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and p~i(t) = 0. The main result is the expression of the solution as the distribution of a composed 
process represented by the standard, homogeneous Poisson process N(i), t > with a random time 
argument T 2u {t),t > as follows: 

K(t)=N(T 2v (t)), t>0. 

The process T 2v (t\t > (independent of N) possesses a well-known density, which coincides with 
the solution to a fractional diffusion equation of order 2v (see (2.15) below). In the particular case 
where v = 1/2 this equation coincides with the heat-equation and the process representing time is 
the reflected Brownian motion. 

These results are reconsidered here, in the next section, from a different point of view, which is 
based on the use of the Generalized Mittag-LefHer (GML) function. The latter is defined as 



oc 

E 

r=0 



ill z r 
r!r(ar + /5)' 



a,/3,7€C, Re(a),Re(f3),Re(j) > 0, 



(1.6) 



where (7) = 7(7 + 1)...(7 + r — 1) (for r = 1,2, and 7 7^ 0) is the Pochammer symbol and 
(7)0 = 1. The GML function has been extensively studied by Saxena et al. (in [25]- [26]- [27] 
and [28]) and applied in connection with some fractional diffusion equations, whose solutions are 
expressed as infinite sums of (1.6 1. For some properties of (1.6 1, see also [29]. We note that formula 
(|1.6[) reduces to ( 1 . 1 ) for 7 = 1 . 



By using the function (1.6) it is possible to write down in a more compact form the solution to 



(1.5 1, as well as the density of the waiting-time of the fc-th event of the fractional Poisson process. 



As a consequence some interesting relationships holding between the Mittag-LefHer function ( 1.1 1 



and the GML function (1.6 1 are obtained 



Moreover the use of GML functions allows us to derive an explicit expression for the solution of 
the more complicated recursive differential equation, where two fractional derivatives appear: 



'PA 



dt 2 » 



2A 



d v p k 
dt v 



-\ 2 (p k -Pk-i), fc>o, 



(1.7) 



for v £ (0,1]. As we will see in section 3, even in this case we can define a process governed by (1.7 1, 
which turns out to be a renewal. The density of the interarrival times are no-longer expressed by 
standard Mittag-LefHer functions as in the first case, but the use of GML functions is required and 
the same is true for the fc-th event waiting-time. 

An interesting relationship between the two models analyzed here can be established by observing 



that the waiting-time of the fc-th event of the process governed by ( 1.7 1 coincides in distribution with 



the waiting time of the (2fc)-th event for the first model. This suggests to interpret our second model 
as a fractional Poisson process of the first type, which jumps upward at even-order events A 2 k and 
the probability of the successive odd-indexed events A 2 k+i is added to that of A 2 k- Correspondingly, 
the distribution of this second process N v (t),t > 0, can be expressed, in terms of the processes N 
and T 2u , as follows: 



Pr 



{#„(<) = fc} = Vr{N{T 2v {t)) = 2k} + Pr{N(T 2u (t)) = 2fc + 1}, k > 0. 



We also study the probability generating functions of the two models, which are themselves solu- 
tions to fractional equations; in particular in the second case an interesting link with the fractional 
telegraph-type equation is explored. 



3 



For v = 1, equation (1.7) takes the following form 



^+2A^ = ~A 2 ( Pfc - Pfc _ 1 ), fc>0 (1.8) 

and the related process can be regarded as a standard Poisson process with Gamma-distributed 
interarrival times (with parameters A, 2). This is tantamount to attribute the probability of odd- 
order values A2k+i of a standard Poisson process to the events labelled by 2k. Moreover it should 
be stressed that, in this special case, the equation satisfied by the probability generating function 
G(u,t), t > 0, < u < 1, i.e. 

d 2 G(u,t) n BG(u,t) , o, . „, „ „ „. 

d \ 2 ' +2A K m ' = X 2 (u-l)G(u,t), 0<^<1 (1.9) 

coincides with that of the damped oscillations. 

All the previous results are further generalized to the case n > 2 in the concluding remarks: the 
structure of the process governed by the equation 

A JIu^T, +■■■+[ _ , A -— = -A (pfe-Pfe-i), A;>0, 



(1.10) 

where ^ G (0,1], is exactly the same as before and all the previous considerations can be easily 
extended. 



2 First-type fractional recursive differential equation 
2.1 The solution 

We begin by considering the following fractional recursive differential equation 



d"p k 
dt v 



-A(pfe-pfe-i), k > 0, 



with p-i(t) — 0, subject to the initial conditions 

P*(Q) = 



k = o 
fe > l 



(2.1) 



(2.2) 



We apply in (2.1l the definition of fractional derivative in the sense of Caputo, that is, for m G N, 

(* 

for v = m 



L„fY> — J r(m-i/) Jo (t-s)i+"-"-> d s ™ u ( s )d s 
17 I £»(t), 



for m — 1 < v < m 



(2.3) 



We note that, for v = 1, (2.1l coincides with the equation governing the homogeneous Poisson 
process with intensity A > 0, so that our first result generalizes the well-known distribution holding 



in the standard case, i.e. Pk(t) 



_ (At)* -At 



k > Q,t> 0. 
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We will obtain the solution to (2.1 1-( 2.2 1 in terms of GML functions (defined in ( 1.6 1) and show 
that it represents a true probability distribution of a process, which we will denote by J\f v (t), t > : 
therefore we will write 

p»(t)=PT{Af u (t) = k}, k>0A>0. (2.4) 



Theorem 2.1 The solution p%{t), for k = 0, 1, ... and t > 0, of the Cauchy problem (|27T])-([2T2[) is 
given by 

(2.5) 



pl(t) = (Af) k E k v tl +1 (-W), fc>0,<>0. 



Proof We take the Laplace transform of equation (2.1l together with the condition (2.2 1 and 
consider that, for the Laplace transform of the Caputo derivative, the following expression holds: 



C 



~df 



u(t);s\ = 



d v 

e- st ^—u{t)dt 



(2.6) 



t=0 



where m = \ v\ + 1. Since, in this case, f € (0, 1] , we get m = 1. Therefore we get the following 
recursive formula, for k > 1, 



while, for fc = 0, we obtain 



since equation (2.1 1 reduces to 



with initial condition po(0) = 1. 



C{ P » (t);s} 

d»Po = 
dt v 



By applying (2.7) iterative ly we get 

C{pl(t)-s} = 



s u + \ 
-Ap ) 



A^s"- 1 



(2.7) 
(2.8) 
(2.9) 



(a» + X) k+1 

which can be inverted by using formula (2.5) of [22], i.e. 

C{t^El 1 {^)-s} = w — ? ., 



(2.10) 



(2.11) 



(where Re(/3) > 0, Re(j ) > , Re(S) > and s > \w\ ) for f3 = is, 5 = k + 1 and 7 = z/fc + 1. 
Therefore the inverse of (2.101 coincides with (2.51. □ 



Remark 2.1 As a cross check, we note that, for v = 1, formula {2.5) reduces to the distribution of 
the homogeneous Poisson process since At) = e~ At jk\. 
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For any v G (0,1], it can be easily seen that it coincides with the result, obtained by a different 
approach in [2], 



by noting that 



r=0 



r\T(isr + vk + 1) 
k\ ^ r!r (v(k + r) + IV ' 

r— 

r\ {-iy- k {\t v y 



r—k 



(r + k)\ = {k+ 1 +r- l)(fc + 1 + r - 2) • • • (k + l)k\ 
= (k+l) r k\. 



Remark 2.2 The result of Theorem 2.1 shows that the first model proposed by Mainardi et al. [14-] 
as a fractional version of the Poisson process (called renewal process of Mittag-Leffier type) coincides 



is expressed in terms of successive derivatives of Mittag-Leffier functions as 



with the solution of equation (2.1) and therefore with \2.5ty . In the paper cited above the distribution 

ives of Mittag-Leffier functions as 

kE u ,i(x) , k>0,t>0 (2.13) 



«*(*) = 



~kT 



dx 



and is obtained by means of the fractional generalization of the Erlang density 



dx k 



E v< i (x) 



(2.14) 



which represents the distribution of the waiting-time of the k-th event. Clearly, for v = 1, fk(t) — 
is the Erlang density representing the waiting-time of the k-th event for the standard Poisson 



(k-iy. 

process (for A = 1). Since we can rewrite (2.13) as 

«*(*) 



Jd 



E 



dx k Z_. T{yj + 1} 



k\ 



E 

j=k 



x — — t u 

jU-i)...(j-k + i)(-t-y- k 



T(vj + 1) 



t*" ~ (i + k)\(-t v y 



E 



A:! ^ l\T(vl + vk + 1) 



it is evident that it coincides with (2.5), for A = 1. 



G 



We derive now an interesting relationship between the GML function in (2.5) and the solution 
to a fractional-diffusion equation, which is expressed in terms of the Wright function 



W atf} (x) = J2 



k=0 



kW(ak + (3) ' 



a > -1, (3 > 0, x e 



Let us denote by v 2v = V2 V (y,t) the solution to the Cauchy problem 



d 2v v _ \2 d 2 v 
dt 2 " ~ A dy 2 > 

v(y,0) = 6(y), 
v t (y,0) = 0, 



t>0, y e R 
for < v < 1 
for 1/2 < v < 1 



(2.15) 



then it is well-known (see [17], p.142) that the solution of (2.151 can be written as 

1 



V2 V {y,t) 



2\t v 



W-uA-u 



m_ 
\t v 



If we fold the above solution and define 

V2v{y,t) 

then, for k > 0, we get 

Plit) - 



0, y<0 



t>0, ye 



y>0 



+OC 



e v -jjv 2 v(y,t)dy 
Yv{N{T 2v {t)) = k). 



(2.16) 



(2.17) 



(2.18) 



In (2.181 T2 V (t),t > represents a random time with transition density given in ( 2. 16 1-( 2. 17 1 and 
independent of the Poisson process N(t),t > (note that from now on N denotes a standard Poisson 
process with intensity A = 1). This result is analogous to what happens for other kinds of fractional 
equations, such as the diffusion ones (see [18]), the telegraph-type fractional equations (in [17]) and 
even the higher-order heat-type equations with fractional time-derivative (see [1]). 

Formula ( 2.18[ ) was obtained for the first time in [2] and we prove it here in an alternative 
form, by resorting to the Laplace transform. We compare (2.101 with the Laplace transform of the 
distribution Pk(t), t > of a standard Poisson process with intensity A > 0, which reads 

X k 



£{ Pk (t);s} = 



( S + A) 



fc+i 



Formula (2.101 can be consequently written as 

C{pl{t);s} = s»- 1 L{p k {t) ] s v }, 



which, by inversion, leads to the following convolution: 



Plit) 



r(i - v) Jo 



(t-w)- u £- x 



- sV y/ x p k (y)dy;w\dw, 



(2.19) 



(2.20) 



(2.21) 
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where p k (t),t > represents the distribution of the Poisson process with intensity A = 1 . 
The inverse transform in (2.21 1 can be expressed as follows 



Cr 



-f-oo 



e s " v/X Pk{y) d V,w \ = / g v {w\ -)p k (y)dy, 



-\- DO 



by recalling that 



Jy/X = I e- sz g v {z; |)cte, < v < 1, w > 0, 



(2.22) 



(2.23) 



where g v (-\y) is a stable law S v (n, (3,o~) of order v, with parameters ji = 0, /3 = 1 and cr = 



Jcos^ 



By inserting p.22\ into < \2.21\ we get 
p£(*) ' 



r(i - y) J 

-\-oo 



(t-w)- 



-Hoc- 



l 



r(i - v) Jo 



y 

g v {w, -)p k (y)dy ) dw 



(t-w) v g u (w;^)dw)p k (y)dy. 



(2.24) 



We recognize in (2.24) the fractional integral of order v of the stable law g v and we apply the result 
(3.5) of [17], which can be rewritten, in this c ase 3 as 

— y 

— . , (t-w) y g v (w;-)dw = v 2v (y,t), 

r(i - v) j a 



where V2v(y,t) is given in (2.171. 
We can conclude from (2.5 1 that 



E k v +i +1 (-\n 



k\X h +H v i h + 1 ) Jo 



- y y k W^ v (-^)dy. 



(2.25) 



Remark 2.3 Formula (2.18) can be useful also in checking that the sum of p k (t) (either in the form 
(2.5) or (2.12)), for k > 0, is equal to one, since it is 

oo „ +oc / OO fc\ „ +00 

$>*(*)=/ e_3/ ETrr 2i/(y '* )dy= / «2,(y,*)% = i. 

fc=0 ^° \fc=0 ' / ^° 

Moreover, since result (2.25) holds for any t > 0, we can choose t = 1, so that we get, by a change 
of variable, 

<tU(- A ) = fc[ I e- x yy k W^ v (-y)dy. 

This shows that the G ML function E„W +1 can be interpreted as the Laplace transform of the function 

iW- v ,x- v (-y). 

Ln particular, for v = ^, since (2.18) reduces to 
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pr{M /a (*) = *} - I ^^f d y ( 2 - 26 ) 

- Pr{JV(|B A (*)|) = *}, 
where B\(t) is a Brownian motion with variance 2X 2 t (independent of N), we get (for t=l) 



i r+°° P -y 2 /4 



The previous relationship can be checked directly, as follows 

(-A) r 2 fe+r r fr + k : 1 



j OO 



v/^fc! ^ r! V 2 2 

= [by the duplication formula] 
1 f (~A) r T (r + fc) (r+fc) 

If R) r (r+fc)! k+1 
k\^ Q H r(r±* + l) h 

2.2 Properties of the corresponding process 

From the previous results we can conclude that the GML function E*+\ +1 (-\t u ), k > 0, suitably 
normalized by the factor (\t") k , represents a proper probability distribution and we can indicate it 
as Pr{jV„(i) = fc}. 

Moreover by (2.181 we can consider the process J\f v (t),t > as a time-changed Poisson process. 
It is well-known (see [8] and [11]) that, for a homogeneous Poisson process N subject to a random 
time change (by the random function A((0, t])), the following equality in distribution holds: 

JV(A((0,t])) = M(t), (2.28) 

where M(i), t > is a Cox process directed by A. In our case the random measure A((0, f]) possesses 
distribution v% v given in (2.16 |-( |2.17 l and we can conclude that J\T V is a Cox process. This conclusion 
will be confirmed by the analysis of the factorial moments. 

Moreover, as remarked in [2] and [15], the fractional Poisson process Af v (t),t > represents 
a renewal process with interarrival times Uj distributed according to the following density, for 
3 = 1,2,...: 

fi(t) = Pr{Uj e dt} jdt = Xt u - 1 E VtU (-\t 1 '), (2.29) 

with Laplace transform 

£{/{•(*);*} = sl 7^ r (2.30) 
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Therefore the density of the waiting time of the k-th event, T k — J2j=i Uj, possesses the following 
Laplace transform 

= (2-3i) 



Its inverse can be obtained by applying again (2.111 for /3 = v, 7 = vk and uj = — A and can be 
expressed, as for the probability distribution, in terms of a GML function as 

f k {t) = Pr{T fe e dt}/dt = \ k f k - l E k VjVk {-\t v ). 



(2.32) 



Formula (2.321 coincides with (2.14 1, for A = 1. The corresponding distribution function can be 

(2.33) 



obtained in two different ways. The first one is based on (2.32) and yields 
F%(t) = Pr {T k <t} 
= X k 



t 00 

vk 



3=0 



j\(k- l)W(vj + vk) 



~ (fc- i + j)i(-\t»y 



v j^j\{k-l)\{k+j)T{vj + vk) 



(k-i + j)\(-xt-y 

^ Q j\(k-l)T(vj+vk + l) 



= \ k r k E* vk+1 (-\i»). 



We can check that (2.33 1 satisfies the following relationship 

PT{T k <t}-Pr{T k+1 <t}=pZ(t). 



(2.34) 



Indeed from (2.33) we can rewrite (2.34) as 



x k t vk E k !Uk+1 (-xtn - x k+ H^ k +^E k +i k+l)+1 (-xtn 



j\(k - l)\T(uj + vk + 1) jlkW(vj + vk+v+l) 

= [by putting I = j + 1 in the second sum] 



{k + m-xn 3 



X k t uk j2 



(k-i+m-xry , M ^ {k+i-m-xn 1 



3=0 
00 



1 x k.vk \ "* \^ f 

j\(k - l)!r(i/j + ^fc + 1) ^ (i - l)!fc!r(i/i + j/Jfe + 1) 



= A V fc 2 - 



(fc + j)!(-A^)J 
j!fc!r(i/j + vk + 1) 



: Pfc(*)- 



The second method of evaluating the distribution function resorts to the probabilities given in 
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the form (2.121: 



Px{T k <t} = Y,p v m {t) 



Finally if we rewrite (2.331 as 



m—k 

oo oo 

EE 

m=k r=m 



r \ (-iy- m (\t v y 

mj r {yr + 1) 



E 



{-\) r (\t v y 



. T{ur+l) 

r—k m—k 



— Y 
i) z - 



(-1)' 



Afcrfc gp-A(-i) fe (-^P 



j=k 



T{vj + 1) 



(2.35) 



and compare it with (2.351, we extract the following useful combinatorial relationship: 



3-1 
k- 



(-l) m , J > k. 



Remark 2.4 As pointed out in [5] and [23], the density of the interarrival times (2.29) possess the 
following asymptotic behavior, for t — > oo: 



Pr {Uj G dt} /dt 



\t v - l E u ^{-\n = - -E vA (-\t v ) 



dt 



(2.36) 



= A 



r iv + 2 r v cos{vn) + 1 
sin (i^r) T(v + 1) v 

A^ +1 = Ar(l - v)V+^ 



dr 



where the well-known expansion of the Mittag-Leffler function 

sin(i/7r) <- + °° 



E v>1 {-\t v ) 



r v-l e -\ 1,v rt 



o r 2jy + 2r jy cos(i/7r) + 1 



dr 



(2.37) 



has been applied (see the Appendix for a proof of (2.37\ l). The density ( 2.36) is characterized by fat 
tails (with polynomial, instead of exponential, decay) and, as a consequence, the mean waiting time 
is infinite. 

For t — > the density of the interarrival times displays the following behavior: 



Pr {Uj G dt} /dt 



Xt 



v-l 



T{v) 



(2.38) 



which means that Uj takes small values with large probability. Therefore, by considering (2.36) and 
(2.38) together, we can draw the conclusion that the behavior of the density of the interarrival times 
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Jers from standard Poisson in that the intermediate values are assumed with smaller probability 
than in the exponential case. 



Remark 2.5 The distribution junction (2.33) of the waiting-time of the k-th event coincides, for 
A = 1, with the so-called Mittag-Leffler distribution of [13] -[19] (see formula (1) of [13]) 

T(t + j)x u W> 



3=0 



jW(t)T(l + v(t + j))- 



x > 0, 



(2.39) 



when the time argument t takes integer values k > 1. On the contrary the space argument x 
coincides, in our case, with the time t. Therefore the process X u (t) with distribution (2.39) can be 
considered as the continuous-time analogue of the process representing the instant of the k-th event 
of the fractional Poisson process j\f v . 



Remark 2.6 We observe that also for the waiting-time density {2.32) we can find a link with the 
solution to the fractional diffusion equation (2.15). This can be shown by rewriting its Laplace 
transform (2.31) as 



£{fZ(t);s} 



A 



£{f k {t);s v } 



(k-iy. 



e~ xt dt. 



By using again ( '2. 23\ we get 



+oo 



9v{t;y)fk{y)dy 
^ (t; A } 0^1)!^ 



(2.40) 



Formula (2.4-0) permits us to conclude that fj((t) can be interpreted as a stable law S v with a random 
scale parameter possessing an Erlang distribution. 

2.3 The probability generating function 

We consider now the equation governing the probability generating function, defined, for any < 
u < 1, as 



(2.41) 



fc=0 



From (|2.1|) it is straightforward that it coincides with the solution to the fractional differential 

d v G(u,t) 



equation 



dt v 



X(u-l)G(u,t), 0<v<l 



(2.42) 
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subject to the initial condition G(u, 0) = 1. As already proved in [2] the Laplace transform of 
G v = G v {u,t) is given by 



£{G„(u,i);s} = / , 
Jo 



f G„(u,t)<it 



(2.43) 



s"-X(u-l) 

so that the probability generating function can be expressed as 

G v (u,t) =E v>1 {\{u-l)t v ). 



(2.44) 



By considering (2.44) together with the previous results we get the following relationship, valid for 
the infinite sum of GML functions: 

oo 

^(^n'KU+ii-^) = e v ,i(Hu- i)o. 



(2.45) 



fe=0 



For u = 1 it shows again that J^fcLnPfeW = 1- The resu lt (2.451 can be checked by resorting to the 
Laplace transforms and noting that 



k=0 



k=0 
c.u-1 



(s v + X) k+1 



(2.46) 



s u - \(u - 1) ' 



Formula ( 2.45 ) suggests a useful general relationship between the infinite sum of GML functions and 
the standard Mittag-Leffler function: 

oo 

X^WJi+iC-*) = - !)): < U < 1. 



(2.47) 



k=0 



For u = 1 it shows again that SfcloPfcW = ^■ 



By considering the derivatives of the probability generating function (2.441 we can easily derive 
the factorial moments of N u (t) which read 



E [M v {t){M v {t) - \)...{N v {t) -r + l)] = 



(Ar) r r! 
T(ur + 1)' 



(2.48) 



These are particularly useful in checking that Af„ (t),t > represents a Cox process with directing 
measure A. Indeed, as pointed out in [11], the factorial moments of a Cox process coincide with the 
ordinary moments of its directing measure. We show that this holds for N v , by using the contour 
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integral representation of the inverse of Gamma function, 



E[A((0,i])] r = 



y r V2 V {y,t)dy 



II 



+00 r 

V 

Xt" 



W- v ,l-u 



y 

Xt v 



dy 



+ OC 



1 1 

X r 

2™ J Ha 
\ rxvr 



z l+vr 



y r dy 



dz 



Ha 

+00 



t »r w r e -w dw 

ATr! 

-dz = 



Ha z^ r(!/r + l)' 



which coincides with (2.481. 



By applying formula (2.47) it is easy to obtain also the moments generating function of the 
distribution, defined as 



fc=0 



which is given by 



(2.49) 



(2.50) 



k=0 



E^(X(e-» - l)t"). 



The r-th moments of the distribution can be obtained by successively deriving (2.491 as 

d r 



EAf v (t) r = (-1) 
= (-1) 



E 



dfi r 
d r 
d[i r 

c jt r (wy 
\ I>J + 1) ' 



/x=0 

E„,i(A(e-"-l)f) 



/i=0 



where Cj jr are constants. For r=lwe get 



m v {t) 



= Xt" 



KAf„(t)=-Xr 
d 



dfi 



dx 



E v>1 (x) 



E V!l (X(e-»-l)t v ) 



Xt v 



(2.51) 



11=0 



x=0 i> + i) ! 



which coincides with the renewal function evaluated in [15], for A = 1. It is evident also from (2.51 1 
that the mean waiting time (which coincides with lim^oo t/m u (t)) is infinite, since v < 1. 



14 



3 Second-type fractional recursive differential equation 



3.1 The solution 



In this section we generalize the results obtained so far to a fractional recursive differential equation 
containing two time-fractional derivatives. We show that some properties of the first model of 
fractional Poisson process are still valid: the solutions represent, for k > 0, a proper probability 
distribution and the corresponding process is again a renewal process. Moreover the density of the 
interarrival times display the same asymptotic behavior of the previous model. 
We consider the following recursive differential equation 



d* v p k 



2A 



d v Pk 



dt 2v ' " dt 
where v € (0, 1] , subject to the initial conditions 



A 2 Ofc-Pfc-i), k > 0, 



Pfc(0) = 
P*(0) = 



k = 
k > 1 



for < v < 1 



(3.1) 



(3.2) 



0. 



k > 0, for - < v < 1 
- 2 



and p-i(t) = 0. In the following theorem we derive the solution to (3.1 1- ( 3.2 1 , which can be still 
expressed in terms of GML functions. 



Theorem 3.1 The solution p\{t), for k = 0, 1, ... and t > 0, of the Cauchy problem (3.1)-(3.2) is 

given by 



'M = >? k t 2kv El%l +1 {-\t») + \ 2k+ h^ k+1 >El k ( f k+1)v+1 {-\n, k>0,t>0. (3.3) 

Proof Following the lines of the proof of Theorem 2.1, we take the Laplace transform of equation 
(3.1 1 together with the conditions (3.2 1, thus obtaining the following recursive formula, for k > 1 

A 2 



C{pr k (t);s} = 



s 2v + 2\s v + A 
A 2 



(3.4) 



(a" + A)- 



£{plU(t);4 



while, for k = 0, we get 



£{%(*); S } = 



s 2 " + 2As y + A 2 ' 
Therefore the Laplace transform of the solution reads 

rt ^ u . , A 2fc S 2 ^ 1 + 2A 2fe + 1 ^- 1 



(s» + A) 



2fe+2 



(3.5) 



(3.6) 



We can invert (3.6) by using (2.11 1 with S = 2/c + 2, /3 = and 7 = 2fci^ + 1 or 7 = (2fc + \)v + 1, 

(3.7) 



thus obtaining the following expression 
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We prove now the following general formula holding for a sum of GML functions: 



x n E^ nv+z {-x) +x n+1 E™ {n+1)v+z (-x) = x n E™- v \ z (-x), n,m>0,z>0,x>0, (3.8) 
which can be proved by rewriting the l.h.s. as follows: 



(to — 


1)! 


x n 




(to — 


1)! 


x n 




(to — 


1)! 


x n 




(to — 


2)! 



E 

3=0 
oo 

E 

3=0 



(m-l+j)l(- X y x r ' 



E 

3=0 



E 

i=i 

oo 

E 



(m - 1 + j)\(-x) j 
j\F(yj + nv + z) 

(m+l-2)\(-x) 1 



x 



(to - 1)! (I - l)\T(vl + nv + z) 

(m + l - 2)\(-x) 1 x 



E 

m-l + l 
I 



(to - 1 + j)\(-x) j+1 



(m+l-2)\(-x) 1 



x" 

T{nv + z) 



= x n E 



m— 1 



(-X) 



For m = 2k + 2, z = l,ar = Ai" and n — 2k formula ( |3.8[ ) gives the following identity: 



which coincides with the first term in (3.3). 

It remains to check only that the initial conditions in (3.2 1 hold: the first one is clearly satisfied 
since it is, for k = 0, 



m) = E 



^E 



(r + l)(-A) r ^ r+1 ) 



-gr(i/r+l) ^ I>r + z, + l) 
= [for f = 0] = 1 



and, for k > 1, 



A 2fc (2fc + r)!(-A) r t"( 2fc+r > A 2fc+i ( 2fc + r + i)!(-A) r ^( 2fc + r + 1 ) 

Pfe ^ ~~ J2k)l r\T(vr + 2kv + l) + (2fe + l)! Hr(z/r + 2fci/ + v + 1) ' 



which vanishes for t — 0. The second condition in (3.2) is immediately verified for k > 1, since it is 



d _A»_ ^ (2fc + r)!(-A)''W 2fc+r )- 1 A 2fc+1 ^ (2fe + r + 1 ) ! (_ A )r r (2fe+r+i)-i 



(2k)\ r\T(vr + 2kv) (2fc+l)! 



r=0 



rT(vr + 2kv + v) 



(3.9) 

which for t — vanishes in the interval | < f < 1. Then we check that this happens also for k = 0: 
indeed in this case (|3.9h reduces to 



- n,'/-i 1 A E 

= 1 v ' r=0 

(_ X y t vr-l Xt u-1 

T{vr) T{v) 



(r + l) 2 (-A) r ^( r + 1 )- 1 



E 

r=2 



T(yr + f) 

(r + l) 2 (-A)T( r + 1 )- 1 ( A^" 1 

r=l 



^E 



T(yr + v) 



I» 



= [for t = 0} = 0. 



1G 



□ 



Remark 3.1 The solution (3.3) can be expressed in terms of the solution (2.5) of the first model as 
follows 

m)=P2 k (t)+P2 k+ i(t)- (3.io) 

Therefore it can be interpreted, for k — 0, 1, 2, as the probability distribution Pr |tV„(£) = fc| for 
a process A/"„. Indeed, by (3.10), we get 



so that it is 



{#,,(*) = fc} = Pr {N u {t) = 2k} + Pr {N v {t) = 2k + 1} 

OO 

5>r{#,(t)=*} 

fe=0 

oo oo 

= Yl Pr ww = 2fc > + E Pr ww = 2fc + 1} = i 



(3.11) 



(3.12) 



fc=0 



fc=0 



Moreover the relationship (3.10) shows that the process governed by the second-type equation can 
be seen as a first-type fractional process, which jumps upward at even-order events A2k while the 
probability of the successive odd-indexed events Aik+\ is added to that of A^k- 



We can check that expression (3.101 is the solution to equation (3.1 , subject to the initial 
conditions (3.2 1, by using the form of p v k appearing in the last line of (2.121 which is more suitable 
to this aim: 



Pk 



;w = E 



r=2k 



2kjT(vr+l) 



E 

r=2k+\ 



{-\ry 



i \2k + l)T{vr + 1)' 



(3.13) 



The fractional derivatives of (|3.13|) can be evaluated by applying the definition (|2.3|), as follows: 
d" 



dt v 



Kit) 



— y 

r(i - v) ^ 



r=2k 



r \ (— A)*Vr 



2k)T{vr + l) J (t-s) v 



-ds - 



(3.14) 



E 



r \ {-X) r vr /"* s^- 1 



ds 



(—\) r t lJr ~ v 



1 ( r\ T{vr + l)T(l - v) 

r(l - v) ^ 



E 



2kJ T{vr - v + l)r (yr + 1) 

r(i/r + l)r(l - i/) 



r(i - 1/) r _^ +1 v 2fc + V r (^ r - v + !) r + !) 



-(-A) r t I/r - I/ 
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and, for v £ [0, i] and m = 1, 



-)2// 



%(*) 



- 2'^ E 
1 



r \ (-\) r vr /"* s"" 1 



T(l - 2v) ^ \2k /T(isr+ 1) J (t - s) 2 » 



ds 



(3.15) 



r(i-2z/) 



E 



1 



r=2fc+l 
oo 



{-\) r vr [* s vr - x 



2k + 1J T(ur+ 1) 7 (i- s) 2y 



-<is 



r\ r(i/r + l)r(l - 2v) 
T(l - 2z/) £^ V 2fc / r(i/r-2z/+l)r(i/r + l) 



r ±vr— 1v 



Y{\-2v) 



E 

r=2fc+l 



r(i/r + i)r(i -2i/) 

2fc + 1/ r(^r - z/+ l)r(z/r + 1) 



(-A)' r ^ 



It is easy to check that the last expression is valid also in the case v £ an d m = 2. By 

considering (3.14) and (3.151 together we get 



d 2v d u 

(-A) r ^ 



(3.16) 



E 

r=2k 



+2A 



2fc/r(w-2i/ + i) 



E 



(-A) r ^ 



E u, 



E 



(-A)''^'-^ 
/2k)Y(vr - v+ 1) 

r + 2\ (-l) r A r+2 r r 



E 



2fc + lyT(^r-2^+l) 



r=2fc-2 



r=2fe+l 

oo 



2k J F(vr + 1) 



E 



r=2k- 



2k+l)T(vr-u+l) 
r + 2\ (-l) r \ r+2 t vr 



\2k+\J T(vr+l) 



-2A 



E 



r=2fe-2 

OO 



E 

L r=2fc-l 

r + 2 
2k 



r+l\ (-l) r X r+1 t vr 
2k J T(vr+l) 



E 

r=2fc 



+ 1 \ (-l) r A r+1 ^ r 



2k + I) T(vr + l) 



E 

=2/c-l 



r + 2 
2fc+ 1 



j 4 



r=2fc-l 



r + 1 



2 e r;>+»E 



= 2A' 



r + 1 
2fc + 1 



^4 , . 



for A r = rTCr+i) — where the second step follows by putting r = r — 2 in the first two sums and 
r' — r — 1 in the second ones. 
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We want to show that (3.161 is equal to 



~>?(Pk -Pk-i) 

(-A)T 



-A 



E 



-2k 



r=2fe-2 
oo 



E 

=2j 

xo 

E 

=2* 
oc 

E 



2fc/ T(vr+ 1) 



(-A) r i" 



E 



r=2fc+l 



2k-2JT{vr+l) 



r=2k 



=2k-2 



A r + 

r 

2k -2 



E 

r=2fc+l 



r 

2fc+ 1 



E 

r=2k-l 
A r - 



{ — \) r t vr 

2k + I J T(vr+ 1) 
(-A) r t" 



2fc - 1/ Y{vr + 1) 



E 

=2fe-l 



r 

2fc- 1 



Then, by considering together (3.161 and (3.17), we get 

-F k + 2\— P »k + x 2 (K-K~i) 



OO 

E 

r=2k-2 
oo 

+ E 

r=2k 
oo 

E 

r=2fc+l 



dt" 
r + 2 
2k 



r + 1 
2k + 1 

r + 2 
2 A; 



r 

2fc- 2 



r 

2k 

r + 1 
2fc 



. 1 , 



2fe-l 



E 

=2 

X) 

r=2fc+l 



r + 2 
2k + 1 



r 

2fc+ 1 



r 

2fc- 1 



r + 2 
2k + 1 



r + 1 
2fc 



r + 1 
2ft + 1 



r 

2fc-l 



r 

2fc-2 



r 

2k +1 



A r + R 



where 



R 



^ 2fc -2 [1-1]+A 2fc _i 

+^2fc 



-A 



2k 



2k + 2 

2k 
2k + 2 
2k + 1 



2fc+ 1 
2fc 

2fc 
2k -2 
2k + 1 
2fc 



-A 



2fc-l 



2fc- 1 
2fc-2 

2fc + l 
2fc + l 



2k 
2k- 1 



,4 



2k 



2k 

2k 
2k + 1 
2fc + l 



2ft- 1 
2fc™ 1 
2k 
2k 



0. 



The sum appearing in (3.181 can be developed by considering that 

(r+1)! 



r + 2 
2k 



r+1 
2k 



(2fc)!(r + 1 - 2ft)! 
r + 1 
2k- 1 



r + 2 



r + 2 - 2ft 



(3.17) 



(3.18) 



(3.19) 
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and analogously 



r 

2k- 1 
r + 2 
2k + 1 



r + 1 
2fc 

r + 1 
2fc + 1 



r 

2k 
r + 1 
2k 



so that, by considering (3.191, we can rewrite (3.181 as follows 



E 

r 

2 k 

E 

r=2fe+l 



r+1 
2fc- 1 



r 

2k + 1 



r 

2fc- 1 



r 

2fc - 2 

r+1 
2fc 



r+1 
2 k 



r+1 
2fc+ 1 



r+1 
2k + 1 



r 

2fc+ 1 



. 1 , 



since 



r + 1 
2Ar — 1 



r 

2k -2 



If we consider now the first two terms of (3.20 1 we get 



so that we get 



E 



since 



r 

2k- 1 



r + 1 
2fc + 1 



r+1 
2fc 



r+1 
2k + 1 



r 

2fc- 1 



r 

2fc+ 1 



A r = 0, 



r 

2fc + 1 



(3.20) 



(3.21) 



3.2 The probability generating function 

As we did for the first model we evaluate the probability generating function and we show that 
it coincides with the solution to a fractional equation which arises in the study of the fractional 
telegraph process (see [17]). 

Theorem 3.2 The probability generating function G„(it,i) = X^fcLo uk Pk(t), < u < 1, coincides 
with the solution to the following fractional differential equation 



d 2v G(u,t) , ox 0"G(u,t) 



dt 2v 



2A- 



dt v 



\ 2 {u- l)G(u,t), 0<v<l 



(3.22) 



subject to the initial condition G(u, 0) = 1 and the additional condition Gt(u, 0) = for 1/2 < v < 1. 
The explicit expression is given by 



G v (u,t) 



lit -X- 1 a / 7/ — 1 

^ -£„,i(-A(l - Vu)t») + ^-^E vA {-\{\ + Ju)t v ). 



2^L 



2^L 



(3.23) 
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Proof By applying the Laplace transform to (3.221, we get 

{s 2 " + 2\s' J )C{G„{u, t);s) + (s 2 ^ 1 + 2A,s ,y - 1 ) = X 2 (u - l)C(G v (u, f); s) 

and then 



C(G v (u,t);s) = 



2Xs h 



s 2v + 2Xs v + A 2 (l - it)' 



(3.24) 



We can recognize in (3.221 the fractional equation satisfied by the characteristic function of the 
fractional telegraph process studied in [17] (see formula (2.3a) with c 2 /3 2 = A 2 (l — u)) and thus the 



Laplace transform (3.24) coincides with formula (2.6) therein. By applying the result of Theorem 
2.1 of the cited paper, we obtain the inverse Laplace transform of (3.241 as given in (3.23). □ 



Remark 3.2 As a first check we note that {3.23) reduces to one for u = 1, so that we prove again 
that (3.12) holds. Moreov er, as an alternative proof of the previous theore m, w e can show that the 
series expansion of (3.23) coincides with 2fcLo u Pk(t) f or Pk(t) given in (3.3): 



G v {u,t) 

Ju~ + 1 (_A(1 - y/u)t")' 



fu- 1 y (-A(l + y/u)t v )i 
2VS Wi + 1) 



(-xty 



Vu + i ^ {xt v y s^[A(/7-,\ k ( 1 \3-k , 



oo 



oo 



2v^ 
AZ+ 1 



fc=0 
oo 



£r E 



u,2mv-\-l 



/u + 1 



s. 2 ^, 1 , , (-xn + E^Li £ (v^o 2m+1 <7 2 + m 2 +1)I/+1 (-A^) + 



vA« - 1 

~2 



m=0 
oo 



E^^(^Ar) 2m ^ +1 (-A^) + 



m=0 



2y^ 

A- 1 

2^/u 



m=0 

oo 



E(-^) 2m+1 ^ 2 + ro 2 +1>+1 (-A^) 



m=0 



oo 

E " m [( A ^) 2m ("AO + (Ar) 2m+1 <7 2 t 2 +1)l/+1 (-At w ) 
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3.3 Properties of the corresponding process 

We can prove that N v (t),t > represents a renewal process, by showing that, also for this model, 
the required relationship between p%.(t) and distribution function of the waiting time Tfe of the k-th 
event holds: 

p£(t) = Pr {f fe < tj - Pr < t} (3.25) 

where 

f fc = inf{i>0:7V;(t) = fc}. 



or alternatively for the Laplace transform of (|3.25l 

C 



{p£(t) ; S } = ~c {my, s} - - s c {Tum, s} 



(3.26) 



In view of relationship (3.101, we can infer that each interarrival time Uj is distributed as the sum 
of two independent interarrival times Uj of the first model and therefore from (2.30) we have that 

A 2 



e~ st Pr 



n 



{% g dt} 



e~ st Pr j#j G dt} 

al- 



ts* + \y 



(3.27) 



By recursively applying (3.261, starting with (3.271, we arrive at 

X 2k 



{fk (*);«} 



(s« + A) 2fc ' 



By applying again (2.11 1 for (3 = v, 7 = 2i^/c and cj = — A we invert (3.281 and obtain 



Moreover from (3.271 or (3.291 we easily get 



Pr jWj G dt} /dt = ZT 1 I 



A 2 



(a" + A)" 



(3.28) 



(3.29) 



(3.30) 



Therefore, in this case, both the waiting-time of the fc-th event and the interarrival times possess 
distributions which are expressed in terms of GML functions. 



Remark 3.3 As a check, we derive directly the probability density of the Uj from that of Uj given 
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in (2.2£ty , without the use of Laplace transform as follows: 



K(s)K(t-s)ds = x 2 



l E u , v {-Xs")(t - sY^E^i-Xit - s) v )ds 

1 00 r \^S7 XsV) \( t -sr-^\- s ^ ds 



a 2 EE 

3=0 1=0 



00 00 



o T(vj + v) 

{-xry+ l 



T(vl + v) 



3=0 1=0 



T(vj + vl + 2v) 



^ T{vm + 2v) 



m=0 



= xH^E^ v (-xn = f^t) 

This confirms that one out of two Poisson events are disregarded in this case (as described in Remark 
3.1). 



Remark 3.4 We evaluate now the asymptotic behavior of the interarrival-times density (3.30), as 
follows: 



00 

Pr |z/, £ dt\ /dt = X 2 ^- 1 



3=0 



(-At")j(j + 1)! 
j\T(vj + 2v) 



(3.31) 



i-i 



(=0 



l(-Xt v ) 
T(vl + v) 



Xt v d 
v dt 



E v , v {-Xt v ) 



1 — v d , „, t d 2 

—E v i(—Xt) + -—E^i-Xt"). 

v at v dt z 



By applying (2.36) and (2.31), we finally get 

Pr [Ui edt\/dt = ^zl x i/u^M_ 

X 2 H sin (z/tt) 



+00 



r 2v + 2r ly cos(i/7r) + 1 



dr + 



+00 



r u+l e -X 1/ "rt 



v ir j r 2v + 2r v cos(i^7r) + 1 
sin(^Tr) 2T(v + \) _ 2v 

\t v + 1 ~~ Ar(l - i^+i' 



rir 



(3.32) 



//we compare \ 3. 32) with the analogous result {2.36) obtained for the first model, we can conclude 
that the interarrival-times density displays the same asymptotic behavior, with the power law decay 
(3.32), so that again the mean waiting time is infinite. 

For t—yQ, instead of (2.38), we get in this case 



Pr jz?,- €dt\ /dt 



XH 2 "- 1 

W 
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The behavior near the origin of the density of the interarrival times Uj has a different structure for 
v < 1/2 (tends to infinity) and v G (1/2, 1] (vanishes for t — > + ). For v = | we have instead that 

Pr jWj E dij is constant for t = and eattaZ £o A . 



We can conclude that Af u is a renewal process and the corresponding renewal function is given 



by 



m„(t) = E/V„(t) = A 2 t 2 "£ Wj 3„ + i(-2At B ) 

z-Eu,v+x{— 2At ). 



(3.33) 



2T(i/ + l) 



Formula (3.331 can be obtained either by deriving the probability generating function (3.231 
(for u = 1) or by using the well-known relationship between the Laplace transforms of the renewal 
function and the interarrival-times density: 



C {m(t); s} 



Indeed, in this case, it is 



C{fh v (t)\s} 



1 



1 C{h(t);s} 

s\-c{h(t)- s y 



s i-c 

A 



{/?(*);«} 



1 



s 1 - 



A 2 



2A 



which gives (3.33 1, by applying (2.11 ) for /3 = v, S = 1 and 7 = 2i/ + 1. 



Remark 3.5 comparing the second form of {3.33) with (2.51), we can note that the following 
relationship between the renewal functions of the two models holds: 



This can be alternatively proved by applying (3.10) as follows: 

00 oo oo 

m v (t) = £W) = EW) + E fe ^+i(*) 

k=0 k=0 k=0 

1 oo 1 oo 1 oo 

= - £(2*)j&(t) + - J2( 2k + - 9 E*w*) 

A;=0 Aj=0 
1 oo 1 oo 

= 9 EiPi(*)-oEftVi(*) 



(3.34) 



(3.35) 



fc=0 



fe=0 



2 ~ 2 E^2fc+i(*)- 



fc=0 



24 



The last term in (3.35), which coincides with the sum of the probabilities of an odd number of events 
of the first model, can be evaluated as follows: 

oo oo 

£*W*) = E A2fc+1 ^ (2fc+1) ^^ + i )+ i(-^) 



k=0 



k=0 

= Xt' y E u ^ + i(—2Xt l/ ), 



as can be checked by resorting to the Laplace transform: 

~ A 2fe+l s ^l A 



fe=0 



( 8 u + A )2fc+2 



s(s y + 2A) 
C{\t v E v>u+1 (-2Xt v );s}. 



Formula (3.34) confirms that the mean waiting time is infinite (as we have noticed in Remark 



3.4): indeed it is 



lim — - — = lim — — lim — - — £^„ +1 (-2Ai ) = 0, 

t — >oo £ t — >oo Zt t — >0 ° Z 



where the second limit can be evaluated by the following considerations: 



Xt 



v-l 



-Ev,v+i{~^Xt v ) 



1 

¥ 
1 

It 



1 - E v ^(-2Xt v )] 
7T 2AF 



(see (5.2) in the Appendix). 



3.4 The special case v — 1 



We consider now the previous results in the special case v = 1. Equation (3.1l reduces in this case 
to the second-order equation: 



d 2 Pk . ^dp k 
and the corresponding solution (|3.3[) is given by 



-A 2 (p fc -Pk-i), k>0 



Pk(t) = 



{Xtf\_ xt , (Xtr+\ _ 



xt 



(2k)\ 



(2k+l)\ 



k>0. 



It is easy to check can show that (3.371 solves (3.361 with initial condition 



Pfc(0) = 
ft(0) = 



1 k = 
k> 1 

k > 



(3.36) 



(3.37) 



(3.38) 
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and P-i(t) = 0. Indeed it is 

y2k£2k-l 



dpk 
dt 

dt 2 

so that we get 



(2A-1)! 
(2fc-2)! ( 



-At 



-\l 



\2fe+2^2fe+l 

(2fc + l)! ' 

\2fc+:U2fc-l 

(2fc-l)! ' 



-At 



-At 



\2fc+2^2fc 



-At 



\2fe+3^2/c+l 

(2fc+l)! ' 



-At 



d 2 Pk . ox dp k 

y2.kj.2k-2 ^ y2k+l£2k-l ^ y2k+2£2k ^ 



2k+l£2k-l 



(2k- 2)' 

2X 

+ (2Ar— 1)! 

^2k^2k-2 



(2k - 1) 



(2k)\ 



O \2fe+3j-2fc + l 
-At ZA 1 -At 
e t^t; n — e 



(2fc + l)! 

\2fc+lj.2fc-l \2fc+2j.2fc 
_ P -At , A t p-At _ A _ p ->* 

(2fc-2)! (2fe-l)! (2fc)! 



^2fe+3^2fe+l 

(2fc + l)! ' 



\2fe+3^2fe+l 

(2fc + l)! ' 



-At 



-A/ 



= -A 2 (p fc -pfc-i). 



Formula (3.371 agrees with the relationship ( |3.10 1 given in Remark 3.1, which in this case we 
can write as follows 

Pk(t)=P2k(t)+P2k+i(t). (3.39) 
Therefore it can be interpreted, for k = 0,1,2,..., as the probability distribution Pr |iV(i) = fcj 

for a "second-order process" N(t),t > linked to the standard Poisson process N by the following 
relationship: 

Pr = k\ = Pr {N(t) = 2k} + Pr {N(t) = 2k + 1} . (3.40) 



From (|3.29|) and (|3.30|) we can easily see that for this process the densities of the interarrival times 

(3.41) 



and of the fc-th event waiting time are given respectively by 

&(t) = X 2 te- Xt 

and 



AW = 



\2fei2fe-l 

(2k- I)!' 



,-At 



(3.42) 



Therefore, in this case, the random variable Tj, representing the instant of the j-th event, is dis- 
tributed as Gamma(X,2j). 

We derive equation (3.371 in an alternative way, which is similar to the construction of the 
standard Poisson process, by considering the relationships (3.111 or, equivalently, the property of 

(3.43) 



the interarrival times described in Remark 3.3. We can write that 



A law „ 
2 j - i 2j, 



where 



T 3 =mf{t>0:N(t)=j} 



2G 



which represents the time of the j-th event of N(t). 



Theorem 3.3 The probability distribution of the process N(t),t > described by (3.401 and (3.431 
solves equation (13.36 1. 



Proof Let us consider now the following intervals 



A = p2k_2,22fc_i) 

B = [T2fc-i,T 2 fc) 

C = p2fc,X2fc + i) 

D = p2fc+l, T 2 k+2) 



so that [T 2fe _2, T 2k+2 ) = Au BUCU D (see Fig.l). 



D t- 



1 2k+l i 2k+2 



Figure 1: The interval [T 2k _ 2 , T 2 k+2) 
We evaluate the following probability, by stopping the approximation at the second-order terms: 



p k (t + 2At) 

,P,(i»( t )-*-M^)*<^(l-^ 



Pr (N(t) =k-l,teA\ (XAt - A 2 (At)'' 



+2Pr (N(t) 



k - 1, t e B ) ( AAt - A 2 (At) 2 ) 1 - AAi 



A 2 (Ai) 2 * 



+ 2P r (Ar ( t) = fc - MeB) ^( 1 _AA t+ A2 ^ 



(3.44) 



27 



Pr (N(t) = k - 1, t G B\ (\At - A 2 (At) 
Pr (#(*) = M G C) (l-AAf+ A ' (At)2 



+2 Pr (iV(t) = jfc, t e c) yi- XAt + A2 -y^- j (aA* - A 2 (At) : 

+ Pr (R(t) = k,teD)(l- XAt + A2(At) j + o ((At) 2 ) , 

where we used the following well-known approximations valid for the standard Poisson process: 

A 2 ( At) 2 

Pr (0 Poisson events in At) = e~ AAt ~ 1 - AAt H K -^- + o ((At) 2 ) 

Pr (1 Poisson event in At) = AAt e~ AAt ~ AAt - A 2 (At) 2 + o ((At) 2 ) 



A , A 2 (At) 2 \ a j A 2 (At) 2 //A 

Pr (2 Poisson events in At) = — \_J_ e -^t _ — \_±_ + Q 



2 \ 



By ignoring the terms of order greater than (At) the probability (3.441 can be rewritten as follows 



p k (t + 2At) (3.45) 

A 2 (At) 2 



2Pr \ N(t) = k-l,t G A 
+ Pr (iV(t) = k - 1, t G A) A 2 (At) 2 + 
+2Pr fjV(t) = i-l,teBj ^AAt - 2A 2 (At) 

+2Pr (N(t) = k- l,t G B 



A 2 (At) 2 



Pr (N(t) = k — l,t € b\ A 2 (At) 
Pr (ft(t) = k, t G cVl - 2\At + 2A 2 (At) 



2 I A+\ 2 

T 

2 / a .a 2 



+2 Pr (iV(t) = jfe, t G (7) (AAt - 2A 2 (At) : 
+ Pr (N(t) = k, t G D \ (l + 2A 2 (At) 2 - 2AAt) 2 + o ((At) 2 ) . 



Since we can write 

Pr (iV(t) = fc — 1, i G A \ + Pr(iV(t) = k - l,t G B) = pfc-i(t) 

and, analogously 

Pr (N(t) = k,teCj + Pr(iV(i) = k,t G D) = p fc (f) 
28 



formula (3.451 becomes 



p k (t + 2At) ~ 2p fc _i(t) 



A 2 (At) 2 



+ p fe _!(i)A 2 (At) 2 + 



+ p fc (t)(l-AAt) 2 +p fc (t)A 2 (At) 2 + 

+2Pr (iV(i) = fe — 1, f G b\ (AAt ~ 2A 2 (At) 2 

+2 Pr f = Me^j (AAt - 2A 2 (At) 2 ) + o ((At) 2 ) 
We consider now the following probability, on a single interval of length At 



p k (t + At) ~ Pr ^JV(i) = k - 1, t G A 
Pr (jV(t) =jfe-l,t£5j ^AAf 

Pr (N(t) = k,teC^j ( 1 - 



A 2 (At)' 



A 2 (At) 2 ^ 



A 2 (At) 2 



+ Pr I N(t) = k,t £ D) I 1 — AAt 



X 2 At 2 



-Pr (Jv(t) = fc-i,tes) 



A 2 (At) 
2 



((At) 2 ) 



(3.46) 



(3.47) 



Pr 



Pr 



(iV(t) = fc-M6flj |^ AAt 

(iv(t) = k,t ec) n - 



A 2 (At)' 



A 2 (At) 2 



+ 



Pfc(t)-Pr(JV(t) = Me<7 



Pfc-l(t) 



1 - AAt 



A 2 At 2 



'((At) 2 ) 



A 2 (At)' 



+ Pr (N(t) =k-l,teB) (AAt - A 2 (At) 
A 2 At 2 \ 



+ p k (t) M - AAt + — ^— J + Pr (%) = Mec) (AAt - A 2 (At) 2 ) + o ((At) 2 ) 



We multiply (3.471 by 2(1 — AAt) and ignore the terms of order greater than (At) , so that we get 



Pfc-i(t)A 2 (At) 2 ~ 2(1 - XAt)p k (t + At) - 2Pr fjV(t) = k- l,t G El) (xAt - 2A 2 (At) 2 ^ 



2p*(t) 



(1 - AAt)' 



A 2 At 2 



2Pr (N(t) = k, t G C) (AAt - 2A 2 (At) 2 ) + o ((At) 2 ) , 
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which can be substituted in the first line of (3.461. Finally we obtain 



p k (t + 2Ai) ~ p k ~i(t)\ 2 (At) 2 + 2(1 - XAt)p k (t + At) + 



(3.48) 



2Pr (N(t) 



l,teB) [XAt - 2A 2 (At) 2 ) -2p k (t) 



(i - xAty 



X 2 At 2 



+ 



2 Pr (N(t) = k,teCj (XAt - 2A 2 (At)'' 



2Pr 



(%) 



k-l,t£Bj (xAt-2X 2 (Aty 



+ 2Pr(iV(t) = /s,te C 1 ) (AAt-2A 2 (At) 2 

+ p k (t) (1 - AAt) 2 + p k (t)X 2 (At) 2 + o ((At) 2 ) 
= p k ^(t)X 2 (At) 2 + 2(1 - XAt)p k (t + At) -p k (t) (1 - XAtf 
We divide (|3.48| by (At) 2 , so that we get 



p k (t + 2 At) - 2p k (t + At) + p k (t) 
(At) 2 



■2XAt m+ ^\: m) = -A 2 (At) 



(At) 2 



2 Pk (t) -Pk-l(t) 

(At) 2 ' 



By letting At — > we easily obtain equation (3.361 



□ 



Remark 3.6 This special case is particularly interesting because it describes a generalization of the 
Poisson process which have been used in many articles. Random motions at finite velocities spaced 
by this particular renewal process have been considered by different authors ([5]-[6]-[10]). 

In particular in [3] it has been studied a model with uniformly distributed deviations which take 
place at even-order Poisson events and therefore its interarrival times are distributed as Gamma(X, 2). 



4 Conclusions 

The results of the previous sections can be generalized to the n-th order case, if we consider the 
following equation 



d nv p k fn\ d^-^pk 
dt nv VI/ dt^' 1 )" 



where v € (0, 1) , subject to the initial conditions 



Pk(0) = 



1 k = 







k > 1 ' 



for < v < 1 



(4.2) 



dti 



Pk(t) 



1 



= j = l,...,n — 1, k > 0, for — < v < 1 



of (4.1l 



and P-i(t) = 0. Following the same steps as for the first two models we get the Laplace transform 

£{^(i); S } = A"£{j55U(t); S }, 



s nv + [ J J A S ("-^ + ... + f i _ x ) A™" 1 * 2 " + A" 
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which can be solved, in view of (2.6 1 and the initial conditions (4.2), recursively, yielding 



£{p*(*);*} 



Y, n -i ( n .)s' /j - 1 \( k+1 ') n - j 

(s^ + A)^ 1 )" 



(4.3) 



For n = 2, we obtain from (4.3) formula (3.6 1. The Laplace transform can be inverted by applying 
again (2.11 1 and the solution is given, also in this case, as a sum of GML functions as follows: 



— V 

Pi 



M = E 



(At 



v,vn(k+l)-vj + l\ AL >■ 



(4.4) 



For n = 1, we get the distribution of the first model (2.5 1, while, for n = 2, we get the solution 
of the second-type equation in the form (3.7 1. In this case the use of (3.8 1 requires much harder 
calculations. Nevertheless a relationship similar to (3.101 can be obtained, even for n > 2, by 
studying the density f%(t),t > of the waiting time of the fc-th event T k . As already seen in section 
3, the following identity must be satisfied by f%(t) : 



C 



{%(*); s} = -C {/£(*); s] - \c {r k+1 (t); s} 



so that, by substituting (4.3 1 in the l.h.s. of (4.5 1 we get 

X nk 



{/*(*);«} 



(s v + X) nk ' 

which can be inverted as usual, thus obtaining 

7m = \ nk t nvk - x E^ vk {-xt v ). 



(4.5) 



(4.6) 



(4.7) 



Again the process is a renewal one, since (4.7) coincides with the sum of k independent and identically 
distributed random variables Uj's (representing the interarrival times) with density given by 



PT{u j edt}/dt = C- 1 { {sV ^ x)n -,t} 



t\ = X n t 



n_i.ni/-l rpn 



(4.8) 



Formula (4.8 1 shows that each interarrival time of the n-th order case is distributed as the sum of 
n independent interarrival times of the first model. This suggests that the following relationship 
between the corresponding probability distributions holds: 



PfcW = Pnk(t) + Pnk+l(t) + -+Pnk+n-l(t), n > 2. 



(4.9) 



5 Appendix 



We derive the integral form of the Mittag-Leffler function (used in (2.37 1), and some generalizations. 
We start by showing that, for < v < 1, 



sin (z^tt) 



r 2u + 2r^cos(^7r) + 1 



dr. 



(5.1) 
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By applying the reflection formula of the Gamma function 

T(ym + l)T(-vm) 



sin(— TTvm) 



we can write the Mittag-Leffler function as follows 
1 ^ (~t») m sm(nvm) 



E 

oo 

E 



Y(vm + 1) 
1 (-l)™sin(7ri/m) r +0 ° 



Y{ym + l)r(— urn) 



IT ' T(^to + 1) J 

TCI— 7 



dr / e- y y um dy 



1 ^ (-l) m sin(7rr/m) 

7T ^ T(^TO+ 1) In. 
m=0 V ' J{) 



+00 



e -ry y um dy dr 



+oo ^ 
e -ry 



E 



Y(ym + 1) 



-f-oo 



1 

~2~7ri 
1 

sin(i/7r) 



2i 

+00 



e- r VE Vtl {-y v e- mv )dy)dr 



r u i gjiriy r ^ _|_ g- 



(//•. 



(//• 



(5.2) 



From (5.1 1 we can derive the following approximation, for large t: 

sin (vir) T(u) 
Eu,i(-t ) * — jT- 

We present now an expression of the Mittag-Leffler function which permits us to interpret it as 
the mean of a Cauchy r.v. We can rewrite (5.1 1 as 

Eu,i(-n 



i 



sin (i/7r) 



n Jo [r v + cos(^7r)] 2 + sin (vit) 



r"- L e- rt dr 



(5.3) 



1 



+ OC 



sin (vir) 



Jo [r + cos(i'7r)] 2 + sin \v-k) 
E x l-e 



l dr 



-tx 1 '- 



[0,oo ) 



where X is distributed as a Cauchy with parameters — cos(7w) and sin(^7r) . Formula (5.3 1 permits 
us to study the particular case where v = 1, since we can write (5.3 1, by means of the characteristic 
function of a Cauchy r.v., as follows: 



EuA-n 



i 

v— >1~ 



+oo 



+ OO 



e~ rt 5(r - l)dr = e"* = J5i,i(-t). 
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Analogously we can study the case where v = 1/2: from the first line of (5.4) and by applying 
formula 3.896.4, p.480 of [7], we get 



i r + °° 9 / r +ao 



1 

7T 

2 

TT 

2 

7T 

1 /tt 



e r *cos(/3r) 



J — oo 
+oo \ 



+ 00 







e " 



'cos(Pr)dr I d/3 



+oo 



dy, 



which coincides with the form given in [17], for x = —yi. 

We generalize formula (5.1) to the case of a two-parameters Mittag-LefHer function: for < v < 1 
and 0</3<^ + l,we have that 



xl-/3 /-+OC 



sin (f7r) 



[r' y + cos(^7r)] + sin (wk) 



sin (/?7r) ... ... . „ . 

f - cos(^tt) + cos(/3tt) 

sin (i/7r) 



dr. 
(5.5) 



We derive (5.5 1 by starting from the series representation of the Mittag-LefHer function: 



y, (-l) ro r m sin((Vm + /3)tt) 

rn— 



r(i/m + /3) 7r sin((Vm + /3)7r) 



i » ( „ 1)m ^» sin ( ( ^ w+/3 )^) 

_ / w? 7>\ T(l — vm — p)T(ym + p) 

tt ^ T(vm + P) tt v ' v y 

m— 

V — 1 ^ r- sin((i/m + /?)tt) / e- rt r- vm - l3 dr \ e^y^+^dy 

,1—0 °° Liim /-+oo / /-+OC \ 

E rv xm sin((^m + /3)tt) / e - rt / e^y^+^dy dr 

t.1—8 p-\-oo ( /*+oo °° / i \m n Mm ai-nvm+i-nB ^ — iTtvm—i-KQ 



1 /. E ^Lwy ir *J* 
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+1-/3 f°° 



2m 

,1-/3 /-oo 



o Jo 



27Ti J 
±1-0 roo 



r v _i_ f>n^v 



dr 



7T 

til 
7T 



-rt v-p rV sin(7r/3) + sin(7r(/? - ;/)) 



r 2 V + 2 r » cos(7ri/) + 1 
sin(7w) 
(r 1 ' + cos(7rz/)) 2 + sin 2 (7r^) 



dr 



v sin(7r/3) sin(7r/3) cos(7ri/) — sin(7rz/) cos(7r/3) 



sin(7w) 



sin(7w) 



From the previous expression formula (5.5) easily follows and, for = 1, it reduces to ( |5.1[ ). For 
(5 = v we obtain from (5.5 1 



( 1 — 1/ /-+oo 



sin (i/7r) 



[r" + cos(^7r)] + sin (un) 



-dr. 



(5.6) 



As a check of (5.6 1 we can study the limit for v — > 1 : 

il — 1/ /•+<» 



sin (f7r) 



1 — 1/ /•+oo 



{r + cos(vn)] +sin 2 (^7r) 



-dr 



2irv J 

/•+QO 



r " e 
re" 



+ OC 



DC 
+ 00 



e -ir/3-|/3| sin(7rw)-i|8cos(7ri/)^ j ^ r 



re 



-rt 



5(r - l)dr = e - * = E hl (-t). 



For large t, the following approximations follow from (5.5) and (5.6 1: 



sin (i/7r) , z> = /3 

We study now a similar expansion for the Wright function, valid for any < v < 1 and 0^1. 
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By following the same steps as in the proof of (5.5), we have that 



(-ty 



E 



m=0 

,1-0 poo 



m\T(i>m + 0) 



2m 
2m 



dy 



i) 



E 

rn—O 



(-e™ v ) m T(vm + (3) 
i\T(vm + j3) r P+ vm 



dr - e- nf) 



E 

m=0 



(-e- <w ) m r(z/m + /3) 
n!r(^?7i + /3) r' 3+t/m 



,1-/3 <-oo 



2™ 



i 1 -^ f 00 e 



r /3 
ri 



-27T/3- 



r 



/3 



dr 



t l-0 f oo e - rt _ 

— 3-e sin i p 



2i 
sin(7w) 



dr. 
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